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SUMMARY 


The results of numerical experiments are presented in which a posteriori 
estimators of error in strain energy were examined on the basis of a typical 
problem in linear elastic fracture mechanics. Two estimators were found to 
give close upper and lower bounds for the strain energy error. The potential 
significance of this is that the same estimators may provide a suitable basis 
for adaptive redistribution of the degrees of freedom in finite element models. 


INTRODUCTION 


One of the most important problems remaining in finite element research is 
the development of adaptive finite element software systems, i.e., finite 
element computer programs which have a local error estimation capability and a 
capability to increase the number of degrees of freedom selectively such that 
the quality of approximation is nearly uniform over the entire solution domain 
and the error does not exceed some pre-specif ied tolerance. 

Research concerned with the development of adaptive finite element soft- 
ware systems has been underway at Washington University for several years. 

This work has resulted in the development of an approach for improving the 
quality of approximation without mesh refinement. In this approach the number 
of degrees of freedom is increased by increasing the polynomial orders (p) or 
introducing non-polynomial basis functions over fixed finite element mesh divi- 
sions. This process of reducing the error of approximation through the addi- 
tion of new basis functions is called "p-convergence" to distinguish it from 
the conventional approach (called "h-convergence") in which the size of finite 
elements (h) is reduced while the number and type of basis functions for each 
element are fixed. 

The efficiency of p-convergent adaptive procedures has been established in 
a series of numerical experiments, reported in references 1 and 2 and it has 
been shown that the rate of p-convergence cannot be slower than the rate of 
h-convergence (I. Babuska, private communication). In fact, for the vast 
majority of practical problems, p-convergence is substantially faster. 

This and other computational advantages suggest that adaptive finite element 
software systems should be based on p-convergence. 

The key problem is to find suitable estimators of error which would indi- 
cate when and where the number of degrees of freedom should be increased over 
the solution domain. Some estimators have been proposed already: 
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1. Babuska and Rheinboldt developed a local, asymptotic, a posteriori error 
estimator for h-convergent approximations (ref. 3). The estimator requires ele- 
ment level computations only and measures the error in strain energy associated 
with an element. An optimal distribution of the degrees of freedom is obtained 
when this error measure is the same for all elements . 

2. Melosh and Marcal (ref. 4) proposed to measure the specific energy differ- 
ence (SED) , defined as the largest difference over the element domain between 
the computed strain energy density function and the same function evaluated at 
the origin of the elemental coordinate system. This measures the effect of 
modes higher than those associated with the (generalized) constant strain 
states on the distribution of strain energy density within finite elements. 

The criterion for mesh refinement is that SED be approximately the same for 
each finite element. In practical computations SED is approximated by the larg- 
est of the strain energy density differences evaluated at quadrature points only. 

3. Peano et al. (ref. 5) proposed a criterion for p-convergent approximations. 
This criterion is based on the rate of change of the total potential energy 
with respect to higher order displacement modes, evaluated before the higher 
displacement modes are actually introduced. When the rate of change of the 
potential energy exceeds a prescribed tolerance, the stiffness terms corre- 
sponding to the higher modes are assembled and the new system of equations is 
solved. The hierarchic structure of the elemental stiffness matrices permits 
efficient use of block relaxation procedures in obtaining improved solutions. 

In this paper we examine the numerical characteristics of a criterion 
similar to that proposed by Babuska and Rheinboldt in ref. 3, but modified for 
p-convergent approximations. Our study is preliminary in nature and is re- 
stricted to one specific problem. 


ERROR MEASURES 


We have examined two measures for the error in strain energy on the basis 
of a problem in two-dimensional elasticity containing a geometric singularity. 
This problem is typical for a large class of problems in linear elastic frac- 
ture mechanics (fig. 1). The error measures were as follows. First we define 
the ith component of the residual vector, which represents the unbalanced body 
force, as 


r . = G u. . . + (A + G) u. . . + X. 
1 1.33 3 .3i 1 


( 1 ) 


in which 

G and X are Lame's constants; 

u. is the ith component of the displacement vector 

J computed by the finite element method. 

The subscripts range over 1,2. 

X. is the ith component of the body force vector. 


44 



-3- 


One of the measures, to be called the "r-estimator," is defined for the kth 
finite element as 


V a) 




+ r„ ) dA 


( 2 ) 


in which 


R k 

*k 


is a constant, to be determined by numerical experiments; 
is the polynomial order of the displacement approximation 
over the kth, element; 
is the area 'of the kth element. 


The other measure, to be called the "t-estimator", is defined over interelement 
boundaries and external boundaries on which tractions are specified, as the 
square of the unbalanced tractions. Specifically, at the boundary of two ele- 
ments, the vector of unbalanced tractions is: 


t ± (s) = [o^ ^ (s) - a ± j ^ (s) ] n^ (3) 

in which 

(a) 

a. . is the finite element approximation to the stress 
1 '' tensor for the ath element; 
n. is the unit normal to the interelement boundary; 

s is the variable along the element boundary. 

At external boundaries the unbalanced traction vector is the difference 
between the computed traction vector and the applied traction vector. When 
displacement vector components are specified, the corresponding unbalanced 
traction vector component is zero. 

The t-estimator is defined as 


T k <6) 



ds 


(4) 


in which 

g is a constant, to be determined by numerical experiments; 

p^ is the polynomial order of the displacement approximation 

at the kth inter element or external boundary segment; 

T denotes the kth elemental boundary, 
tc 
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Both the r and t measures were found to give close indications of the 
total error in strain energy in p-convergent approximations. The details are 
as follows. 

The sample problem represented in fig. 1 does not have a known exact solu- 
tion. For this reason it was necessary first to estimate the exact value of 
the total strain energy U. This was possible by utilizing the asymptotic 
relationship given in reference 6: 


U = U 


P 



(5) 


in which 


U is the computed (total) strain energy, based on pth 

^ order polynomial approximation; 

c is a constant; 

NDF is the (net) number of degrees of freedom. 


2 2 

Extrapolating on the basis of eq. (5), U was found to be 0.7702 O £ /E 
(+ 0.0002 a 2 & 2 / E) in computations involving two different finite element mesh 
divisions. The computed dimensionless strain energy values, the corresponding 
errors and the values of the two error estimators for parameter values a = 3, 
g = 2 are given in Table I. The percent changes in the estimators as p is in- 
creased tend to bound the corresponding percent change in the strain energy 

error with increasing precision such that the change in the r-estimator is 
smaller and the change in the t-estimator is larger than that of the strain energy 
error. Furthermore, the values of these percentage changes are monotonically 
decreasing for p > 2. This suggests the possibility that two constants c and 
c. could be found such that a relationship. 


c E R, (3) < U - U < c Z T (2) (6) 

r k k " P - ‘ k k 

remains valid for all p-values and that the upper and lower bounds become pro- 
gressively closer to the strain energy error as p is increased. This is, of 
course, highly speculative at the present because no theoretical justification 
exists, but consistent with the observations of this numerical experiment. For 
example, if we choose c = 0.7008 and c = 0.2657, the two estimators will give 
the value of the strain r energy error at p = 7. The resulting relationship 
between the strain energy error and the two estimators is illustrated in 
fig. 2. 

The question naturally arises whether the same estimators would bound the 
energy error at the element level as well. Clearly, this approach will be use- 
ful only if the indicators tend to zero with increasing p at about the same 
rate as the error in strain energy does, not only for the entire solution 
domain but for individual finite elements as well. The presently available 
information is sufficient only to indicate trends. 
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The problem chosen for study does not have a known exact solution; thus 
the energy error cannot be computed with precision. For the global solution 
the rate of convergence formula (eq. 5) provided a basis for extrapolation to 
the limit value of the computed strain energy values. Such formulas are not 
available for predicting convergence at the element level. For this reason an 
ad hoc procedure for extrapolation had to be devised. It was assumed that the 
strain energy converges at the element level as 

u (k) « u (k) + a- (7) 

■ p 

(k) 

where y and v are constants, u is the strain energy of the kth element, and 

u p (k) is the strain energy of the kth element computed on the basis of uniform 
pth order polynomial approximation over all of the finite elements. Taking two 
consequtive values of polynomial orders q and p, y can be eliminated and the 
extrapolated value of u^ k ' is 


v (k) v (k) 

,, , p u — q u 

t (k) = B_ 9 

v v 

p - q 


( 8 ) 


(k) 

where v and u were chosen such that for q = 5, p = 6 and q = 6, p = 7 the 
value of uOO was constant. The resulting estimate of the strain energy for 
element 2 is u( k ) = 0.0844 When the p-distribution is uniform over 

the entire mesh, the estimators vary over the finite elements by several orders 
of magnitude. As could be expected, their value is the greatest for the crack 
tip elements (element numbers 1,3,4) and least for the elements remote from the 
crack tip (element numbers 5,6,8). In those elements which are not on the 
crack tip, the estimators apparently approach zero faster than the error in 
strain energy. 

When the distribution of p is altered such that the element to element 
variation in the estimators is reduced, then the estimators tend to become 
closer to the energy error. Specifically, letting p = 7 for the crack tip 
elements and p = 3 for the remote elements, the variation in the estimators is 
reduced somewhat but is still 5 orders of magnitude for the r-estimators and 4 
orders of magnitude for the t-estimators. In this case we have complete third 
order polynomial approximation in the transition elements 2 and 7 with some ad- 
ditional shape functions ranging in order from 4 to 7. The estimators for ele- 
ment 2 were computed as c r R.2(3) = 0.2 x 10 -5 , c t T 2 ( 2 ) = 0.2 x 10 4 , and the 
dimensionless strain energy error as 0.4 x 10 -4 . Thus the indications are 
that the bounding property of the estimators may be preserved at the element 
level only if the p-distribution is such that the estimators do not vary by 
more than 2 or 3 orders of magnitude. 

In all computational experiments conducted to date the r- and t-estimators 
were found to be consistent indicators of the source of energy error, in this 
case the crack tip singularity. 
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CONCLUDING REMARKS 


A preliminary numerical investigation has indicated that readily computable 
bounds may exist for the error in strain energy in p-corivergent finite element 
analysis. These bounds would provide an indication of where the degrees of 
freedom should be increased over the solution domain in adaptive finite element 
analysis to achieve uniform quality of approximation. 
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TABLE I 

ESTIMATORS FOR THE ERROR IN STRAIN ENERGY 


[Centrally cracked panel, 8-element mesh, uniform p-distribution. *] 


P 


d p U - U P l V 3 > 


S T (2) 
k R 


2 

0.6936 

0.0766 

0.0864 

0.4662 

3 

0.7305 

0.0397 

0.0534 

0.1989 

4 

0.7461 

0.0241 

0.0321 

0.1042 

5 

0.7538 

0.0164 

0.0218 

0.0658 

6 

0.7584 

0.0118 

0.0161 

0.0452 

7 

0.7613 

0.0089 

0.0127 

0.0335 


2 2 

*To convert entries into dimensioned values, multiply by a SL /E. 
^Poisson's ratio: 0.3. 




A C B 


Fig. 1. Centrally cracked square panel mesh division 
and element numbering. 



Fig. 2. Variation of the error indicators and the error 
in strain energy with the polynomial order (p) . 
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